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The distinctive features of the electronic structure of vortex states in superconducting graphene 
are studied within the Bogolubov-de Gennes theory applied to excitations near the Dirac point. We 
suggest a scenario describing the subgap spectrum transformation which occurs with a change in the 
doping level. For an arbitrary vorticity and doping level we investigate the problem of existence of 
zero energy modes. The crossover to a Caroli - de Gennes - Matricon type of spectrum is studied. 
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Recent exciting developments in transport experiments 
on graphene [l[ have stimulated theoretical studies of 
ossible superconductivity phenomena in this material 
0, 0, HI ■ A number of unusual features of supercon- 
ducting state have been predicted, which are closely re- 
lated to the Dirac-like spectrum of normal state excita- 
tions. In particular, the unconventional normal electron 
dispersion has been shown to result in a nontrivial modifi- 
cation of Andreev reflection and Andreev bound states 
in Josephson junctions Q. In this Letter we consider 
another generic problem illustrating the new physics of 
Andreev scattering processes in graphene, namely, the 
electronic spectrum in the core of a vortex that can pre- 
sumably appear in the superconducting graphene in pres- 
ence of magnetic field. The electronic vortex structure 
for Dirac fermions has been previously studied in parti- 
cle physics for a situation equivalent to the zero doping 
limit in graphene, when the Dirac point lies exactly at 
the Fermi level. A number of important results have been 
obtained, namely, the exact solutions for the zero energy 
modes [1, 0] and the subgap spectrum for some model 
gap profiles within the vortex core (lpj . The problem of 
zero energy modes has been further addressed in [ll|, [l2| 
for vortices in various condensate phases described by the 
Dirac theory on a honeycomb lattice. Our goal is to de- 
velop a theoretical description of the electronic structure 
of multiply quantized vortices beyond the zero doping 
limit and study the transformation of the spectrum un- 
der the shift [i of the Fermi level from the Dirac point. 

The energy spectrum is determined by the vortex 
winding number (vorticity) n and labeled by the angu- 
lar momentum v which is a conserved quantity for an 
axisymmetric vortex. For zero doping the sub-gap spec- 
trum consists of n zero energy states [1, [t| . The states 
with higher energies lie close to the gap edge ±|Ao|. In 
the present Letter we demonstrate that, with an increase 
in doping level /_t, the distance between the energy lev- 
els decreases, so that more and more states fill the sub- 
gap region. The set of low-energy levels gradually trans- 
forms into a set of n "spectrum branches" E^ (v) where 
i = 1, ... ,n. If v is considered as a continuous param- 



eter, these n energy branches cross the Fermi level [13[ 
as functions of v. The detailed behavior of the energy 
spectrum as a function of /i crucially depends on the 
parity of the winding number. For odd n there exists 
one branch which intersects zero energy at v = 0. Its 
crossing point belongs to the spectrum, thus resulting in 
an exact zero energy mode. The crossing points of other 
77—1 branches do not generally belong to the spectrum 
for finite [i, thus these zero modes do not exist. However, 
some of these n — 1 zero energy modes can appear again 
for certain doping levels. For high doping ^S> |Ao|, 
they appear almost periodically with an increase in the 
Fermi momentum hkp = \fj,\/vp by a characteristic value 
8k p ~ Here £ is a superconducting coherence length 
£ = Hvf/Aq, and Ao is a homogeneous gap value far from 
the vortex center. In a singly quantized vortex the en- 
ergy spectrum for high doping is given by the Caroli-de 
Gennes-Matricon (CdGM) expression similar to that for 



usual s-wave superconductors 14|. For an even winding 
number the crossing points of all n branches do not gen- 
erally belong to the spectrum if fj, ^ 0. In this case the 
exact zero modes are absent. Nevertheless, some of these 
zero energy modes can appear again for certain fj, in the 
way similar to the odd vorticity case. For high doping, 
they appear periodically with increasing fi. Indeed, the 
energy branches cross the Fermi level at certain momenta 
i/W ~ The zero modes appear when these be- 

come integer (odd n) or half-integer (even n) for discrete 
equally spaced (due to equidistant energy levels) values 
of the Fermi momentum. Because of the symmetry with 
respect to the sign of energy (see below) the zero energy 
modes appear and disappear in pairs, for v and —v. 

Model. The Bogolubov-de Gennes (BdG) equations 
in superconducting graphene can be written as two de- 
coupled sets of four equations each Q . Assuming valley 
degeneracy we can consider only one of these sets known 
as "Dirac-BdG equations" $,UJ^ 
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-A) u + At) = (E + n)u , 
-A J v + A*u = (E - fj,)v . 



(1) 
(2) 
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Here p = — iTiV, and u — (tti, u 2 ) and u = (ui, V2) 
are two-component wave functions of electrons and holes 
in different valleys in the Brillouin zone [lq . The in- 
dices 1 and 2 denote two sublattices of the honeycomb 
structure. They form spinors in the pseudo-spin space, in 
which the Pauli matrices cr Xl a y , a z arc defined; wc also 
introduce a vector & — (a x , a y ). The energy is mea- 
sured from the Fermi level. For zero doping, the Fermi 
level lies exactly at the Dirac point, p = 0. For nonzero 
doping, the Fermi level is shifted by p > upwards (elec- 
tron doping) or downwards p < (hole doping) from the 
Dirac point. For a homogeneous gap A = Ao and zero 
magnetic field the wave functions take the form of plane 
waves it, v oc e ip r and we obtain two noninteracting en- 
ergy branches E 2 — A§ + (p^fVFP) 2 , which correspond to 
the pseudospin orientation chosen parallel and antiparal- 
lel to the momentum direction, respectively. 

Vortex states. Consider an n-quantum vortex, A = 
I A(p)|e m< ^, where p, <f> are cylindrical coordinates. For an 
axisymmctric magnetic field the eigenstates are labelled 
by a discrete angular momentum quantum number v: 



e ffl */ 2 (7(p) 



(3) 



The extra factor e~ l<T ^ 2 as compared to the usual BdG 
functions u and v comes from the angular dependence 
of the momentum operator in cylindrical coordinates 
p x ±ip y = e ±lc ^h[—id/dp±p~ 1 d/d(j)]. The transformation 
properties of the wave functions under rotation around 
the vortex axis correspond to those for an s-wave super- 
conductor with the replacement n — *■ n — 1. As a result, 
v is half-integer for even vorticity n and integer for odd 
n in contrast to standard s-wave superconductors. 

General properties of zero energy states. Equations 
([1]) and ([2]) are invariant under the transformation 



E — » — E , ia y u* 



V , lUyV 



(4) 



Thus, for arbitrary vorticity and doping level, a set of 
zero energy modes (v,i,Vi) labelled by an index i, sat- 
isfies [|| Vi = ia y u* , Hi = —i&yV*. This transforma- 
tion couples the states with opposite angular momenta. 
One can separate two types of zero energy solutions: (i) 
The eigenfunction components transform into each other, 
i = j, which can be realized only for v = 0. (ii) There 
exists a pair of eigenfunctions i =/= j with opposite v 7^ 
coupled by the above transformation. 

Consider solutions of the first type which can exist only 
for an odd-quantum vortex. We put v = ia y u* and find: 



Vfct ■ 
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u + iA& y ii* = 



(5) 



We look for the solution in the form u = CG )^ ; where 

fj(0) 

is a normal state solution of Eq. ([5]) with A = and 
£ is a real function satisfying 



a x e 



'+U<®?to F (d(/dp) = a y U<®*AC 



In a homogeneous magnetic field H the functions 
diverge at large distances p of the order of the mag- 
netic length Lh — \JhcjeH except for a discrete set of 
p which correspond to the Landau energy levels. How- 
ever, similarly to the case of usual superconductors, this 
large-distance divergence is cut-off cither at the magnetic 
field screening length (isolated vortex) or at the intervor- 
tex distance ~ Lh (flux lattice). Thus, we can disre- 
gard the above mentioned divergence for the analysis of 
the states localized near the vortex core at distances of 
the order of the coherence length £. Considering weak 
magnetic fields near the vortex core H <C H C 2 such that 
Lh 3> £ = hvp/A we can even neglect the vector po- 
tential at all. In this way we obtain for p > 0: 



L> (0) - C 



VWV("- 1 W 2 J (n _ 1)/2 (fe F p) > 
e ^/4 e «(»+ 1 )*/ 2 J (n+1)/2 (feFp) 



C(p) = e~ K ^ ; K(p) 



hv 



F Jo 



\A(p')\dp' . (6) 



To address the problem of the zero energy solutions 
of the second type for v 7^ we note that the sym- 
metry transformation Eq. |4]) implies that under the 
transformation E — > — E and v — > — v the functions U 
and V in Eq. (J3J) change according to V2 — ► U\ and 
V\ — * — U-x. Therefore, for E = 0, the functions with op- 
posite momenta coupled by Eq. ((U obey V^ v = U\- v , 
Viu = —U2 —v Equations for U\ and U2 take the form: 
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Ui, v . (8) 



Here N± = u+ (n± l)/2 and AAp) is the <jy component 
of the vector potential. We multiply the first equation 
with U 2 ,v and the second one with Ui jV . Next we do 
the same for the opposite sign of v and then add all the 
resulting equations together. We thus find: 



{Kv - U l-u - Kv + K-u) P d P 



>2,v 

d 



hv F I \p(U ljV U ilU - U lt - v U2,- v )] dp . (9) 
Jo dp 

For a solution regular at the origin and decaying at p — > 
00 the r.h.s. of Eq. © is zero and, thus, the l.h.s. of this 
equation also should vanish. One can see that for v 7^ 
and p ^ 0, the l.h.s. is nonzero in general. Therefore, 
zero-energy levels do not generally exist. 

To illustrate this we analyze the case of small p using 
the perturbation theory. For p = the equation for 
Ui decouples from that for {/ 2 . As follows from [§] for 
positive circulation n > 1, the functions Ux, v are regular 
while U2,v = if — \{n — 1) < v < i(n — 1). For negative 
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circulation, the functions are regular while U\. v = 
if |(n + 1) < v < — i(n + 1). In total there exist |n| zero 
energy levels. For small /i we can calculate the integral in 
the l.h.s. of Eq. © using the wave functions satisfying 
Eqs. © - © with ti = 0. Using the results of Ref. [|| 
for n > 0, one can check that U 2 — U 2 L V is nonzero unless 
v = 0. Therefore, there is no decaying solution and, 
thus, there is no zero energy level for small [i and v ^ 0. 
However, the zero-energy levels may appear occasionally 
for certain finite values of p, as it is shown below. 

Large doping levels. Quasiclassical equations. The 
limit of large /x is very instructive and helps to get the 
complete picture of the spectrum transformation. For 
the analysis we follow a standard quasiclassical scheme 
(see (l8j for details) and introduce a momentum rep- 
resentation: ip(p) — (2nh)~ 2 J rf 2 pe lpp / R V(p); where 
p = p(cos 6> p , sin 6> p ) = pp . Let us transform the 
wave functions so to choose the spin quantization axis 
along the direction of the momentum p: u — Sg u , 
v = Sg v , where S = e~ w " a ^ 2 {a z + a x )/V% and S + = 
(a z + a x )e ie v a * /2 /^/2. Equations H), © take the form: 

[v F a z p - p] g u + h A g u + h A g v = Eg u , (10) 
[p - v F a z p\ g v + h A g v + h\g u = Eg v , (11) 

where we define the operators h A = — (ev F /c)S + &A.S, as 
well as h K = S+AS and h+ = S+A*S. Here A = A(p) 
and A = A(p) are functions of the coordinate operator 
p = ihd/dp. Generally, the operators h A and mix the 
energy bands with opposite pseudospin orientations with 
respect to the momentum direction. However, taking the 
limit of large positive chemical potential p S> A and 
considering only the subgap spectrum one can adopt the 
single band approximation with a fixed pseudospin orien- 
tation: & z g u = g u = g u (l,0), & z g v = g v = g v (l,0). The 
accuracy of such approximation can be determined using 
the second order perturbation theory; the corresponding 
corrections to the energy caused by the off-diagonal pseu- 
dospin terms in h A and areSEA ~ (^ 4 / X| J -)(Ao/ kp£), 
SEa ~ Ao/(fc^^) 3 . These corrections are much smaller 
that the proper energy scale Ao /kp£ for the sub-gap 
spectrum. One concludes therefore that the single-band 
approximation is sufficient for the case of large doping 
p ^> A. The assumption of large /x > allows us to 
consider only the momenta close to hkp. Thus, we put 
p = HIcf + q (\q\ <C hkp) and perform a Fourier transfor- 
mation into the {s,9 p ) representation: 



— kp 



dse- iqs/h g{s, 9p) . (12) 



We introduce vector g and Pauli matrices f x , f y , f z in 
electron-hole space. The angular dependence can be sep- 
arated: g(s,9 p ) = exp(iv6 p + if z n6 p /2)F{s). 

We start with an odd-quantum vortex and consider 
the low-energy levels of the first kind which belong to 



the anomalous energy branch crossing the Fermi level at 
v = 0. For angular momenta vjkp^ <C 1, the solution 
can be found using a linear expansion of h x and in 
terms of the angular momentum operator v — —id/d6 p . 
Assuming a small homogeneous magnetic field, we obtain 



illVFT- 



dF 
~ds~ 



|A( S )| 



F = E'F .(13) 



where E' = E + vhuj c /2 and lo c — eHvp / chkp is the 
cyclotron frequency; v is an integer. Considering the last 
term in the above Hamiltonian as a perturbation we find 



E 



n 
lkp~ 



s 2 



(14) 



where I = / °° ( 2 (s)ds while K(s) and ((s) are defined 
by Eq. ©. Eq. |Q3J) corresponds to the CdGM result 
14j with an account of a magnetic field in the spirit of 



Ref. 17[. Since v can now take zero value, the spectrum 



has a zero energy mode. 

The levels of both first and second kind for any vor- 
ticity n can be easily described in the limit of large v. 
Here we can replace the angular momentum operator by 
a classical continuous variable. Instead of Eq. ([13]) we 
get from Eqs. (|10p . (jlip the Andreev equations along the 
rectilinear quasiclassical trajectories: 



-lhV F T z — 
OS 



hjJrkpb 



g(s,0 p )=O, (15) 



where D = ? x ~ReA(p) — f a ImA(p), b = —vjkp is the 
trajectory impact parameter, the position vector on the 
trajectory is p = (scos6> p — bsm6 p , ssm9 p + bcos9 p ). 
Eq. ((T5|) describes the quasiparticle states in an arbi- 
trary gap profile. For the particular case of a multi- 
quantum vortex the gap operator in Eq. (|15p takes the 
form: A(p) = |A(Vs 2 + b 2 )\[(s + ib)/Vs 2 + b 2 ] n e m9 ". 
For \E\ < Aq the solution of corresponding Andreev 
equations is known to give |n| anomalous energy branches 
E^ (6) crossing zero energy as functions of a continuous 
parameter b (see 13j, [l8( and references therein). 

The total wave function should be single- valued. The 
appropriate Bohr-Sommerfeld quantization rule for the 
angular momentum reads: kp J Q n b(9 p )dO p = 2nm + 
ir(n — 1), where m is an integer. The last term accounts 
for all the phase factors that appear in S and g. As 
a result, the angular momentum v should be integer or 
half-integer for odd and even vorticity, respectively. This 
agrees with the conditions imposed by Eq. ©. The ex- 
pression for the anomalous branches at low energies take 
the form E^{b) - A (6 - &*)/£, where -£ < < f. 
For the branches with bi ^ the change in the Fermi 
momentum is accompanied by the flow of the eigenval- 
ues through the Fermi level: the energy levels cross it by 
pairs at discrete values kp = \u/bi\. 
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FIG. 1: The subgap spectrum vs the angular momentum v 
for a singly (a) and doubly (b) quantized vortex in super- 
conducting graphene with different doping levels: p/Ao = 1 
(circles) and p/Ao = 10 (squares). The quasiclassical CdGM 
solutions of Eq. (|15[l are shown by solid lines (dashed parts of 
lines are guides for eye). We choose here = 100. 

Numerical results. The above analytical considera- 
tions are in excellent agreement with the results of 
our numerical calculations. For the numerical solution, 
Eqs. ((T|), @ are expanded in the basis of eigenfunc- 
tions for a normal graphene disk of radius R. The vor- 
tex is placed at the disk center. We are not interested 
here in the effect of boundaries on the subgap spectrum 
and, thus, for the sake of simplicity we use the condi- 
tions of zero current through the boundary in the form 
ui(R) = vi(R) = 0. To suppress the influence of the 
boundary condition the radius R is taken much larger 
than the eigenfunction decay length in the superconduct- 
ing phase, i.e., the coherence length £. The normal- 
state eigenfunctions satisfying the above boundary con- 
ditions are: U^ p oc Um],p{p), V^ p — and Up = 0, 

and m is either m e or nih, Tn e ^ = v — (1 n)/2, while 
k^R are the p-th zeros of the Bessel function J m (x). In 
a finite-size disk the basis also includes the electronic 
and hole functions: (^ p -Z-i^j for m e> h < — 1, respec- 
tively, which correspond to k e ' h — 0. The gap func- 
tion of n-quantum vortex is approximated by A(p) = 
A e m< ^ [p/ \/p 2 + £ 2 ]™- F° r the numerical diagonalization 
procedure the matrix is truncated keeping the number 
of eigenstates larger than 10i?/£. Shown in Fig. 1 are 
typical energy spectra for singly and doubly quantized 
vortices calculated for different doping levels. Spectrum 
consists of anomalous branches crossing zero of energy 
and branches lying close to the gap edges. Figures 1(a) 
and 1(b) clearly demonstrate the difference between odd 
and even vorticities. Fig. 1(a) illustrates also a crossover 
to a spectrum of the CdGM type with increasing doping 
level. For large p the two anomalous branches in a doubly 
quantized vortex are well described by the quasiclassical 
analytical solutions found in Ref . fl8l] . 



Discussion. We can summarize that in an s-wave su- 
perconducting graphene the subgap spectrum in a vor- 
tex core has a set of energy branches which cross Fermi 
level as functions of v treated as a continuous parameter. 
The number of the branches is determined by the vortex 
winding number n. Whether the crossing points of these 
branches with zero belong to the energy states or not de- 
pends on the parity of the vortex and on the doping level 
p. For a degenerate case p = there are \n\ zero energy 
levels, i.e., all crossing points belong to the spectrum. 
If the doping level p is increased, the spectrum trans- 
formation depends strongly on the parity of the winding 
number, (i) For an odd winding number one of the zero 
energy modes (type I mode), namely that corresponding 
to v = 0, survives with an increase in the doping level. 
The other \n\ — 1 (i.e., even number of) levels split off 
zero with increasing p (type II modes), (ii) For even- 
quantum vortex, all \n\ levels belong to the second type 
and split off zero. In general, the zero energy modes of 
the second type can appear or disappear with the change 
in the Fermi momentum and may exist only in pairs. 
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